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Abstract 

Numerous recent papers have proposed pooled testing strategies to assess the associa- 
tion between a group of rare genetic variants and complex human traits. We consider 
a unified framework that addresses key issues in pooled testing, including the use 
of weights assigned to particular variants and the directions of genetic effects. We 
categorize methods into two classes: linear statistics sensitive to specific directional 
alternatives and "omnibus" quadratic statistics, which have reasonable power across 
a wide range of alternatives. The powers of the statistics are related to non-centrality 
parameters associated with normal approximations. The power for a quadratic statis- 
tic that we consider depends only on the explained variation for a group of variants. 
In contrast, the power for linear statistics depends on the variants' frequencies, di- 
rections of their effects and weights, if used. The complex relationship among these 
factors suggests that, even if rarer variants tend to have larger genetic effects, the 
common strategy that chooses weights inversely proportional to a rare variant's fre- 
quency can adversely affect power. Further, even if the causal variant effects all have 
the same direction, quadratic statistics can outperform linear statistics unless the 
proportion of causal variants in the group is sufficiently high. We also show that 
recent methods that use data-driven weights to maximize linear statistics are opera- 
tionally similar to quadratic statistics. Our framework can deal with categorical and 
quantitative traits, trait-dependent selection of individuals in a study, and it extends 
to include stratification or adjustment for covariates. We discuss the performance of 
tests based on analytical results, simulation studies and applications to GAW17 data. 
We conclude that considerable attention should be given to background information 
that can guide the selection of SNPs for pooled testing. 

Keywords: linear statistics; quadratic statistics; score tests; weighting; power; next 
generation sequencing. 
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1 Introduction 



Genome-wide association studies (GWAS) have identified numerous genetic variants 
(single nucleotide polymorphisms, or SNPs) that are associated with complex diseases 
or traits (e.g. Manolio et al. (2008), Hindorff et al. ( 2009[ )). However, because of their 
limited sample size such studies are effective only at identifying common variants, that 
is, for which the minor allele frequency (MAF) is not too small (e.g. MAF > 5% for 
sample size ~ 2000). In addition, variants that have been identified through GWAS 
explain only small fractions of the genetic components of disease risks or variability 
in traits ( [Manolio and Collins (2009)). There is now mounting evidence that rare 
variants (as represented by SNPs with small MAFs) may contribute significantly to 
phenotype variation but because they are rare, their discovery is more difficult (e.g. 



Li and Leal 


( 


2008 


); 


Bansal et al. 


( 


2010 


); 



scale studies involving huge numbers of individuals might not be a viable option due 
to cost, heterogeneity and other concerns, attention has focused on methods that 
combine information across multiple rare SNPs in a genomic region, based on data 
generated from the next generation sequencing (NGS) technology. This area is the 
focus of our article. 

Papers that propose pooled association testing strategies based on the combination 



of information across multiple SNPs include Morgenthaler and Thilly (2007), Li and 



Leal (2008), Madsen and Browning (2009), Bansal et al. (2010), Han and Pan (2010), 



Hoffmann et al. (2010), Morris and Zeggini (2010), Price et al. (2010), Yi and Zhi 



(2011 


), 


Neale et al. 


(2011) 


Wuet ah 


(2011 


) and 


Lee et al. 


(2012 



provided various solutions but insight into settings when a method will perform well. 



indifferently or poorly is still limited. Recently, Lin and Tang (2011) have given a 



theoretical and empirical evaluation of the methods in the preceding papers, and have 
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proposed a new test statistic based on the maximum absolute value across a set of 



statistics. Lee et al. (2012) have also compared previous tests and give a new family 



which they term SKAT-0, for "sequence kernel association test - optimal". Basu and 



Pan (2011) have conducted an extensive empirical evaluation study of these methods 



in the context of case-control studies. 

Comparisons across many of the proceeding papers do not result in firm conclu- 
sions about which method to use when preparing to address a specific genetic asso- 
ciation testing setting. Unfortunately, the precise prior information essential to an 
informed decision concerning the choice of test statistic is in any case usually unavail- 
able. Moreover, the term "optimal" , used in connection with any of the many classes 
of proposed tests, has meaning in a theoretical framework where precise knowledge 
exists of the true underlying phenomena, but the relevance to real settings is unclear. 



Basu and Pan (2011) reached more specific conclusions, which we discuss later, but 



recommend that when prior information is absent, both linear and quadratic statistics 
(see below) be used. 

In this paper we consider tests for genotype-phenotype association within a uni- 
fied framework. Most test statistics that have been proposed can be divided into 
two classes: linear composite statistics which are powerful against specific association 
alternatives (e.g. Morgenthaler and Thilly (2007), Li and Leal ( |2008 ), Morris and 



Zeggini (2010), Madsen and Browning (2009) and Price et al. (2010)), and quadratic 



statistics that have reasonable power across a wide range of alternatives (e.g. Neale 



et al. (2011), Wu et al. (2011)). A feature of many of the linear statistics and of the 



quadratic statistics of Wu et al. (2011) and Lee et al. (2012) is the use of weights 
associated with individual SNPs, one rationale being that there is evidence suggesting 
that larger effects are associated with rarer SNPs. We study both classes of statistics 
theoretically and empirically, and provide new insights. We consider scenarios that 
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involve varying proportions of beneficial (protective), harmful (deleterious) and neu- 
tral SNPs. We show that although linear statistics can perform well if the majority of 
SNPs are causal and harmful (or causal and beneficial), they perform poorly relative 
to quadratic statistics when there are both protective and deleterious SNPs and more 
generally, where a substantial portion of the SNPs under consideration are neutral. 
We also demonstrate that the effects of using weights in linear or quadratic statistics 
that are inversely proportional to a variant's MAF depend on the type of statistic. 
Even if the assumption that rarer variants tend to have larger genetic effects is true, 
such weights can in some cases have an adverse effect and in others a beneficial effect. 
We do not specifically study statistics using adaptive weight selection; however, we 
show that for linear statistics, adaptive methods of weight selection without external 



information (e.g. Han and Pan (2010), Yi and Zhi (2011), Hoffmann et al. (2010), 



Lin and Tang (2011)) are operationally similar to using quadratic statistics. We also 



consider the question of optimality and indicate why it is in practice unachievable. 
Our discussion deals with all types of traits (categorical, quantitative) and allows 
trait - dependent selection of individuals in a study or non-independent SNPs. It 



also extends to include stratification or adjustment for covariates. Like Basu and Pan 



( |2011[ ), |Lin and Tang| ( [201l| ), |Wu et aL] ( |201lD and |Lee et aL] ( |2012D , we focus on score 
statistics, which are both theoretically and computationally efficient. Our results are 
relatively transparent and easy to apply to practical situations. 

The remainder of the paper is organized as follows. Section 2 introduces notation 
and the framework for testing for association between rare variants and phenotypes. 
We consider arbitrary traits in Section 2 and then the important case of binary traits 
in Section 3. Section 4 presents theoretical power calculations for quantitative traits 
that indicate when various methods will do well, and Section 5 gives simulation results 
comparing different approaches with case-control settings. Section 6 examines data 
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from GAW17 data provided by the 1000 Genomes Project (Almasy et al. (2011), 
1000 Genomes Project Consortium (2010)). Section 7 discusses extensions to deal 
with covariates. Section 8 concludes with some recommendations for pooled testing. 



2 General Traits 

We assume that a group of J SNPs labelled j = 1, . . . , J and a trait Y are under 
consideration. The objective is to consider whether there is association between 
Y and one or more of the SNPs; we do this by formally testing the hypothesis of 
no association. Let Aj and Bj represent the rare and common alleles for SNP j, 
respectively. For a set of n unrelated individuals, let Yi be the measured trait and 
Xij denote the genotype for the i*'^ = 1, . . . , n) person's j*'^ (j = 1, . . . , J) SNP. For 
simplicity we assume that Xij denotes whether the rare allele Aj is present {Xij = 1) 
or absent [Xij = 0) in person i and let Xi = {Xn, Xij)' . It is straightforward 
to consider the case where Xij is the number (0, 1 or 2) of the rare allele for SNP 
j, but since Aj is rare there will be no or very few individuals with genotype AjAj 
in a study of typical size. Tests below can also be readily modified to accommodate 
vectors Xij that code genotypes AjAj, AjBj or BjBj in other ways. We assume for 
now that there is no adjustment for covariates; we address this in Section 7. 

Most proposed methods (e.g. see Basu and Pan ( 2011[ )) for testing a null hy- 



pothesis of no association between Y and X are based on statistics Sj {j = 1, J), 
which individually measure association between Y and given SNP j. Without loss of 
generality we assume that Sj is such that E[Sj] = and Var(S'j) = (Jq^ under the null 
hypothesis; under alternatives, we denote E[Sj] and Var{Sj) by fij and cr|. There 

n 

are many options for 5*^; for example 5*^ = YlO^i ~ is natural choice if the 

i=l 

Yi are approximately normally distributed and is also used with binary traits. The 
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approaches referred to in Section 1 can be expressed in terms of statistics of the form 



^3 ~ y^.^i^ijy 



J=1,...,J 



(2.1) 



i=l 



where ctj is a function of either Yi or its rank, with Y17=i = 0. Such statistics arise 
naturally from regression models relating Y and X, as we discuss later. 
Our interest is in testing the null hypothesis 



Hq : Y and X are independent. 



(2.2) 



We first review the permutation distribution of S = [Si, Sj)'. Although exact or 
asymptotic model-based distributions can be obtained in many cases, the permutation 
distribution is often used to compute p-values; this is the distribution that arises 
from randomly permuting Yi,...,Yn and assigning then to the Xi. Under Hq, the 



permutation mean of S is and the covariance matrix S5 has entries (e.g. Kalbfleisch 



and Prentice (2002), Sec. 7.2) 



T,p{j,e) = covp {Sj, Si 



i=l 



n — 1 
\ / 



n — 1 



mn - 



mjTJii 



n 



(2.3) 



where rrij = J2 Xij and niji = Yl XijXu, for j = 1, . . . , J and £ = 1, . . . , J. This 

i=l i=l 

also applies when F is a discrete variable (see supplementary materials), when the 
genotypes Xij are correlated within individuals (e.g. due to linkage disequilibrium, 
LD) and when sampling of individuals is Y-dependent. 
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Many authors have considered hnear test statistics for Hq of the form 



Wl = ^WjSj = w'S, 



(2.4) 



where the weights Wj are specified positive values and w = {wi, wj)'. Basu and 



Pan (2011) provide a review, but we note two cases: Morgenthaler and Thilly (2007) 



considered the "cohort allehc sums test" (CAST) where each wj = 1, and Madsen and 



Browning (2009) based wj on the population MAF pj, with larger weights for SNPs 
with smaller pj. The rationale for the latter weights is that deleterious SNPs would 
be subject to "purifying selection" and so be rarer in the population than neutral 



SNPs, but evidence for this so far seems slight. Price et al. (2010) also considered 



"threshold" versions in which Wj > only if the relative frequency of SNP j is below 
a specified threshold (e.g. 1% or 5%). Such statistics can have good power against 
alternatives where E[Sj] = fij > 0, with fij > for some subset of {j = 1, . . . , J}. 
However, their power may be poor for alternatives where both positive and negative 
values of fij are possible and, in addition, when only a small proportion of the J 



SNPs have /ij > (Neale et al. (2011), Basu and Pan (2011)). This is studied here 
in Sections 4 and 5. 

To facilitate further discussion, we assume without loss of generality that Y is 
defined so that a SNP with fij > is termed deleterious (harmful) and one with 
Hj < is termed protective (beneficial). A SNP with fij = is termed neutral and 
deleterious or protective SNPs are termed causal. In the absence of prior knowledge 
concerning the effects of SNPs, a preferable family of statistics might be of quadratic 
form Wq = S'AS, where A is a positive definite (or semi-definite) symmetric matrix. 
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A common choice is A = and this is effectively a Hotelhng statistic, 



(2.5) 



where S5 is the covariance matrix of S under Hq. This statistic arises from models for 
Y given X in many settings, and we consider it in theoretical and numerical studies 
below. 



Basu and Pan (2011), Neale et ah (2011), Wu et al. (2011), Lee et ah (2012) and 



others have considered other quadratic statistics. 



Wq = S'AS 



(2.6) 



For example, the 'SSU' statistic (Pan (2009)) and 'C-alpha' statistic of Neale et al. 



(2011) are based on A = I, the J x J identity matrix and the SKAT statistic of Wu 



et al. (2011) uses A = diag{wi, ...,wj), where the wj are weights. Linear statistics Wl 
in (2.4) can also be expressed in this form, since Wl is given by (2.6) with A = w'w. 



Statistics of the form (2.6) can be obtained from random effect regression models 



in which Y is related to X through a linear function f3'X and the J x 1 regression 
coefficient /3 is a random vector with mean and covariance matrix tA. The hypoth- 



esis r = then corresponds to (2.2) and a score statistic for testing it is (Goeman 



et al. (2006), Basu and Pan (2011)) 



Wr 



Q 



-S'AS - -trace(AS< 



Using Wq is equivalent to using (2.6). A number of authors have claimed that as an 



"omnibus" test statistic, the Hotelling statistic (2.5) lacks power in many settings 



Conversely, the fact that (2.6) can be obtained as a score test for single parameter (r) 
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has suggested to many that it will have more power than (2.6). It is often overlooked 



however, that the choice A = S^^ in (2.6) produces (2.5) and so it can be obtained 



as a "single parameter" test. The performance of a test statistic (2.6) depends on A 
and on the mean fj, and covariance matrix of S* under alternative hypotheses and 
we explore this in what follows. 

It is instructive to consider the case where S is normally distributed. The vectors 
S considered here and by others are all at least asymptotically normal, and analytical 



derivations of power and discussions of optimality (e.g. Lin and Tang (2011); Lee 
et al. ( 2012[ )) rely on this. The normal case is well known in multivariate analysis, 
in connection with tests for a multivariate mean /j, = E{S), where fi = (/ii, 



see for example Mardia et al. (1979), Ch. 5. For tests of Hq: fi = 0, the power of 
statistics such as as (2.4)-(2.6) against an alternative hypothesis Hi for which n ^ 
depends on /j, and on the distribution of S under Hi. In particular, suppose that 
under Hi the distribution of S is multivariate normal with mean /j, and covariance 
matrix S. For simplicity we assume S is known; this is the case for some statistics 
and asymptotically it is generally all right to assume it. The following distributional 
results then hold: let Ai, Aj be the eigenvalues of S^/^AS^/^ and let P be the J x J 
orthogonal matrix whose columns are the corresponding eigenvectors. Then 

(i) Wq is distributed as a linear combination of independent non-central random 
variables, 

J 

(2.7) 



where ncj = (P'S ^^^Ai)j and xir denotes a non-central random variable 
with k degrees of freedom and non-centrality parameter r. 

(ii) Under the null hypothesis ^ = 0, Wq is a linear combination of independent 



Xi random variables; each ncj = in (2.7). 
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iii) li A = T, then Wq ~ Xj.nc where nc = ^t'S fi. In this case VFq = Wh of 



(2.5) 



(iv) Zl = {w' Sy / (w'T^w) is distributed as Xi,nc with nc = {w' n)"^ / (w'T^w). 



For such distributional resuhs, see Rao (1973), Sec. 3b. 4. We note that software exists 



for the computation of probabihties associated with hnear combinations of central or 
non-central Xi random variables. In particular, we used the CompQuadForm package 



in R (Duchesne and de Micheaux (2010)). 

The results above allow the power against an alternative hypothesis where S 



N(fj,,T,) to be calculated for any linear test statistic (2.4) or quadratic test statistic 
( |2.6 ). Critical values for a test oi Hq: = are obtained according to (ii). We note 
in particular that 



(a) For a (two - sided) size a test using the linear statistic (2.4) or equivalently Z 



in (iv) above, the a critical value for Zf is ~ '^)) ^^e 1 — a quantile for the 
Xi distribution. The power against an alternative {n, S) is then 



where ucl = {w' / {w'T,w) 



(b) For a size a test using the Hotelling statistic (2.5), the a critical value for Wh 
is ~ ot)- The power against an alternative (^, S) is 



P{x\nc„ > - «)) 



where ncn = A^'S /x. 



A number of authors have claimed to obtain "optimal" tests. This is theoreti- 
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cally possible if we specify a suitable family of test statistics but for this to be of 
practical use we must have strong prior knowledge about the alternative hypothesis. 



For example, among the family of linear statistics (2.4), maximal power is obtained 



when w = H ^/j,; when the Sj are independent so that S = diag{af,...,aj), this 
gives Wj = fij/cr"^ (j = 1, J). This linear statistic is in fact optimal among all tests 



of fixed size based on S. Quadratic statistics (2.5) or (2.6) for which A has rank 2 



or more can never be optimal against a specific alternative (^t, S). However, they 
can maintain high power over wide ranges of alternatives, whereas a linear statistic's 



power can be poor except near a specific alternative. Goeman et al. (2006) and others 



have discussed optimality of score statistics (2.6) coming from random effects models 



but these results are based on averaging over a family of alternatives, which may or 
may not be plausible in a given setting. 



The Hotelling statistic (2.5) is a reasonable choice when alternatives with both 



deleterious {fij > 0) and protective {fij < 0) SNPs are plausible and also performs 
well more generally, as we show later. It can be shown that for a given linear statistic 
Wl = w'S we can decompose Wh as 



Wh = ZUw) + Q{w) 



{21 



where Q{w) and Z^iw) are independent and, under an alternative (a*, S), Q{w) is 
non-central Xj-i ncn with 

1 V 



ncQ = ncH — ncL = /^'S ^fi — {w' / {w'T^w). 



(2.9) 



The linear statistic will have low power when ncq > is large, while it is optimal 
when ncq = 0. 



The test statistic Wc = S'S, given by (2.6) with A = I, has been found powerful 
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in a quite wide range of settings (e.g. Neale et al. (2011), Basu and Pan (2011)). For 
the most part, the settings investigated were ones where the regression coefficients 
I3j in a model for Y given X were unrelated to the frequency (MAF) of the rare 
variant. In cases where causal SNPs (and larger |/3j|) are more likely to be found 
among rarer variants, the situation may, however, be reversed. To illustrate this, 
let pj = P{Xij = 1) in population and suppose the individuals represent a random 
sample (l^,Xj), i = 1 



,n. The covariance matrix in (2.3) is approximately a 



multiple of the diagonal matrix diag{pi{l — pi), ...,pj{l — pj)} when the SNPs are 

J c2 



mutually independent. Then, (2.5 



is approximately Yl 



Pi(i-Pj) 



so that Sj is weighted 



SNP j. It is thus possible that Wh may be more 



inversely according to the MAF 

J 

powerful than Wc = Yl ^f- We investigate this further in Sections 4 and 5. 



Finally, some authors (e.g. Han and Pan (2010), Hoffmann et al. (2010), Lin and 



Tang (2011)) have proposed two-stage or other adaptive approaches in which the 



weight vector w in (2.4) is chosen after preliminary examination of the direction of 



Sj or an estimate of its effect. However, such approach cannot on its own (i.e. without 
the use of background information from other sources) improve globally on quadratic 
statistics. In fact, if we choose the w that maximizes the standardized linear test 
statistic (2.4), then we end up with the quadratic statistic ( 2.5[ ). In particular (see 



e.g. Mardia et al. (1979) p. 127 or Li and Lagakos (2006), Sec. 3) 



'HFXVariWL 



(w'S) 



.-1 



where the maximizing vector is w = ^- This helps explain why 



Basu and Pan 



(2011) found that adaptive procedures did not perform as well as one might hope. 



We remark that Lin and Tang (2011) have proposed a test statistic Tmax based 



on the maximum of a specified set of K linear statistics, each with different weights. 
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T| = {w'j^SY / {w'jjls'^k)- We do not consider such statistics here, but it is clear that 
their performance will depend on the choice of 'appropriate' weight vectors w^. In the 
case where there is little prior information and the Wk are selected to cover a range 
of alternatives with deleterious and protective SNPs, it seems likely that max{T^) 
would be similar to Wh- A similar suggestion involving quadratic statistics is made 



by |Lee et al.| ( |2012[ ). 

In practice there is often very limited prior information about the nature of /-i, 
especially concerning which SNPs might be causal, so one cannot be confident that 



a linear test statistic (2.4) will be effective, nor which quadratic statistics might be 
best . In Sections 4 and 5, we investigate situations in which a specific statistic will 
be more or less powerful. 



3 Dichotomous Traits 

Many previous investigations have focused on binary responses and case-control sam- 
pling. Suppose Yi = 1 and Yi = 0, respectively, indicate whether an individual has 
a certain condition ("case") or not ("control"). A widely used statistic for assessing 



association between binary Yi and Xij is (e.g. Basu and Pan (2011), Neale et al. 



(2011)) 



Sj = ^{Yi- Y) Xij = Tj - niiTij/n, (3.1) 



1=1 

n n n 

where Tj = ^ XijYi, rrij = ^ Xij and ni = ^ Yi. When Xij denotes the presence 

i=l i=l i=l 

{Xij = 1) or absence {Xij = 0) of the rare variant for SNP j, Tj is the number of 
cases with the rare variant and rrij is the total number of individuals with the rare 
variant; rii is the number of cases. The statistic S = (S'l, S'j)' has expectation 
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zero under Hq for either random or case-control sampling; when cases are rare in the 
population, the latter sampling design is typically used. 

Under either type of sampling, the conditional distribution of Sj given mj, n and 
Hi, and under Hq, is hypergeometric, 

where uq = n — rii. This is also the permutation distribution for Tj when the rti cases 
and no controls in the sample are randomly assigned to individuals 1, . . . ,n, and it 
is the basis of Fisher's exact test applied to the 2x2 table defined by Xij = or 1 



and = or 1. In our case we want to consider all J SNPs, and we find from (2.3) 
that the covariance matrix for S has entries 



^p{j,^) = —f TT "^i^ ' 3.3 

n{n — 1) V n 



n 

for j = 1, J; i = 1, J, where ruji = ^ XijX^ is the number of individuals with 



the rare variant for both SNP j and SNP i. The diagonal entries in (3.3) agree with 



variances obtained from (3.2), but the Tj may not be (conditionally) independent for 
j = l,...,J. 

As in Section 2, choice of a statistic for testing Hq should be guided by which 



alternatives for fij (j = 1, . . . , J) are of interest Linear statistics of the form (2.4) 
have been considered in case-control settings by several authors cited above. Recently, 
Nealeetal consider the quadratic test statistic 

W^^t{s^^'^). (3.4) 
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This is, for large n and small mj/n, approximately equal to 



W'i = X^(^|-V^arp(S,)) 

= S'S — trace (Sp) . (3.5) 



This test arises from a random effects model as discussed in Section 2 and is also a 



" C-alpha" (score) test for extra-binomial variation in the Tj. Basu and Pan (2011) 



noted the essential equivalence of Wq to Wq, and to Wc = S'S. 



The distributions of test statistics such as (2.4), (2.5) or (3.4) under null (Hq) 
and alternative (Hi) hypotheses can be approximated in sufficiently large samples by 
normal, chi-squared, or linear combination of chi-square distributions. This follows 
directly from the fact that n~^^'^S is asymptotically normal as n goes to infinity 
with rii/n and rio/n fixed, and an application of results in Section 2. Although the 
approximations may be inadequate in some situations, examination of power under 
(approximate) normality of S provides considerable insights for both quantitative 
and binary traits. Thus we next assume normality for S in Section 4 and provide 
analytical power comparisons for quantitative traits. We then return to case-control 
settings in Section 5 and empirically evaluate the performance of test statistics via 
simulation. 
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4 Power Calculations Under Normality 
4.1 Theoretical calculations 

We consider a setting where Xij denotes whether individuals i has [Xij = 1) or does 
not have {Xij = 0) the rare variant of SNP j, and assume a normal linear model. 

Y^ = ^o + PiXa + ... + PiX.j + a for i = l...n (4.1) 

with Cj ~ A^(0, 0"^) and the Xij mutually independent Bernoulli variables with P{Xij = 



1) = pj {j = 1, J). Explained variation of the J SNPs under model (4.1) i 



IS 



^ Var{E{Y\X)) ^ E/=i - ^ ^ PA^-P,)P] ^ ^ 

when EV is small. We will refer to EVj = pj{l—pj)(3j /a"^ as the "explained variation" 
due to SNP j. The score statistic S = {Si, Sj)' with 

n n 

^, =Y,{Y,- Y)X,^ = J2 - X,)Y„ (4.2) 

i=l 1=1 

arises from maximum likelihood theory for testing Hq : (3 = (/3i, ..-iPj)' = 0. Due 
to normality of the Yi, the distribution of 5*^ given the Xi {i = 1, ...,n) is 

Sj N{Aj/3j,Aja'^) (4.3) 

where Aj = mj{l — rrij/n). For simplicity we consider the case where rrij and mji 
equal their expected values, npj{l —pj) and npjPi {j 7^ /), so that 

5~iV(^,S<,). (4.4) 
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where ^l = npj(l-pj)/3j)' and S5 = diag{npi{l-pi)a'^, ...,npj{l- 

PjW}. 



We consider two linear statistics (2.4) : Wis with weights wj set to 1 and Wlw 



with weights Wj set to 1/ a/ — Pj))- We also consider two quadratic statistics 



(2.6) : Wc with matrix A = I (C-alpha) and Wh with matrix A = E^^ (Hotelling). 



Under model (4.1) we have from the results in Section 2 that 

i) Wls is distributed as Wls ~ NiJ2j=i '^Pji^ ~ Pj)l^j^ Z]/=i " PjO^^^ 

ii) Wlw is distributed as PVny ~ ^(X]/=i ~ Pj)(^jjiT'-J'^'^) 



iii) VFc is distributed as Wc ~ XI -^iXinc 5 where Xj = npj{l — Pj)cr and ncj = 

i=i ' ' 

npj{l-pj)l3]/a^ = nEVj. 

J 

iv) Wh is distributed as W^// ~ Xjnc where nc = ^ npj(l — pj)f3j/a^, which is 
approximately equal to nEV. 

Thus the power of Wh depends (approximately) just on the total explained vari- 
ation and sample size and is not sensitive to the specific effects of the causal variants. 
On the other hand the power of the C-alpha statistic depends on the explained varia- 
tion of the SNPs and also corresponding "weights" Xj. The effect of the Xj is to give 
larger weight to more common SNPs (SNPs with larger MAF). The power of a linear 



statistic Wl (2.4) is a function of the non-centrality parameter 



, . n {Ej=iWjPji^ - Pj)l3j) 
ncL{w) = ^— —— — . 4.5 

^ Ej=iw]Pj{l-Pj) 

To get weights near the optimal Wj = Pj/cr'^ requires significant prior information, 
and the effect of weights inversely proportional to pj is unclear. We provide numerical 
results on power of Wls, Wlw, Wc and Wh under various conditions in the next 
section. 
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4.2 Numerical Comparisons 

In order to consider a broad range of scenarios, we randomly generated 1000 different 
models. Each model has randomly selected values for J, J^, J^, Pj and pj as follows: 
J - number of SNPs was randomly selected from the set {10, 20, 30, 40, 50}; pc — Jc/ J 
- proportion of causal SNPs was generated from U{0.1, 1); p^ = JdI Jc - proportion 
of deleterious SNPs was generated from U (0.75, 1), /3j - genetic effect of SNP was 
generated differently depending on the scenarios described later and pj - probability 
of the rare variant for the SNP was generated from [7(0.005,0.02). Sample size 
was fixed to be n = 1000 and the level of the test was set to be a = 0.0001 to reflect 
the fact that testing would typically be conducted for multiple genetic regions. Since 
the power of W^s-, Wlw: and Wh is a function of values ^j/a and pj, without 
loss of generality we set cr^ = 1. We consider two different situations concerning pj 
and genetic effect ^j. The first scenario assumes no relationship between pj and j3j 
for each SNP and takes \fij\ ~ C/(0.45,0.5) for each causal SNP. Explained variation 
EVj of a single causal SNP ranges between 0.1% and 0.49% under this scenario with 
SNPs with smaller MAP having smaller explained variation. The second scenario 
assumes that explained variation due to a single causal SNP is independent of MAP 
and therefore SNPs with smaller MAFs have larger genetic effects Explained 
variation EVj of a single causal SNP is generated from [7(0.001,0.0025) in this case, 
and the genetic effect is then determined from EVj — ^^^^ ^^^''^^ . 

Figure 1 shows results based on the 1000 randomly generated models under the 
first scenario. It compares the analytically calculated power of linear statistics with 
and without weighting, Wls and Wlw- It also compares power of the quadratic 
Hotelling statistic Wh and C- alpha statistic Wc- In view of the wide variations in 
model parameters, the powers of the tests vary widely across the 1000 models, but that 
the power of the two linear statistics is similar for each model, and the power of the 
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two quadratic statistics is similar for each model. Moreover, neither statistic within 
each class dominates across a majority of the models. Figure 2 has results based on 
the 1000 randomly generated models under the second scenario. We now see that 
for the linear statistics the picture is similar to that in Figure 1, but the Hotelling 
statistic performs better than the C-alpha statistic across almost all models, as our 
earlier comments suggest. 

Figures 3 and 4 compare linear and quadratic statistics. Figure 3 compares power 
of the hnear statistic Wis and the HoteUing statistic Wh for 1000 models generated 
under the first scenario, and Figure 4 compares power of Wls and Wh for 1000 
models generated under the second scenario. In this case we consider three sample 
sizes: n — 500, 1000 and 2000. Figure 3 indicates that with sample size n — 500 
there is low power for both statistics in a large fraction of the 1000 scenarios. The 
linear statistic more often achieves a moderately high power, in essence because only 
models with fairly high proportions of causal SNPs with similar (same direction) 
effects produce much power. As n increases, however, the quadratic statistic displays 
good power across many models and by n = 2000 dominates the linear statistic for 
most of the models. In Figure 4 the linear statistic dominates when n — 500, but 
power exceeds 0.5 in only small proportion of models. When n — 1000 the linear and 
quadratic statistics are best about equally often, but the linear statistic achieves high 
power more often. When n — 2000, however, the quadratic statistic dominates in the 
vast majority of the scenarios. 

We also investigated the relationships between model parameters and the power 
of the linear {Wls) and Hotelling {Wh) statistics. We present results for settings with 
J = 30 and n — 1000 in Figure 5. Results are based on the 10,000 randomly generated 
models under the first scenario. Figure 6 has results based on the second scenario, 
with explained variation of each SNP independent of pj. To show the relationship 
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between power of the tests and the number of causal SNPs we grouped models by the 
number of causal SNPs; to show the relationship between power and the proportion of 
deleterious SNPs, we sub-grouped models by the number of deleterious SNPs. The X 
axis in Figures 5 and 6 indicates both the number Jc of causal SNPs and the number 
Jd of deleterious SNPs in each model. Values Jc — 3, 4, 30 give the main scale in 
each figure, with Jd ranging from 0. 75 Jc up to Jc between successive values for Jc- 
Figures 5 and 6 suggest a number of important conclusions: 

i) Performance of the linear statistic varies widely across scenarios. For Wls to 
perform well it is necessary not only that the effects of causal SNPs are (almost) 
all in the same direction (deleterious or protective) , but also that the proportion 
of causal to neutral SNPs is not too low. 

ii) The quadratic statistic Wr performs well across a range of scenarios with vary- 
ing proportions of deleterious, protective and neutral SNPs. It can outperform 
linear statistics even in cases when causal SNP effects are all in the same direc- 
tion, if the ratio of causal to neutral SNPs is not high. 

iii) The powers of both Wls and Wh strongly depend on the percentage of causal 
SNPs in the group of SNPs. For settings here with n = 1000 and realistic levels 
of explained variation for causal SNPs, powers of 0.5 or more require that a 
majority of the SNPs be causal. The results suggest that considerable attention 
should be given to background information that can guide selection of SNPs 
for pooled testing, and that rather large sample sizes may be needed to provide 
adequate power. 

Finally, we generated Jd/Jc form [7(0.75, 1) here to reflect the common assump- 
tion that rare causal SNPs are more likely to be deleterious than protective. However, 
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we also simulated models where Jd/ Jc was U (0.5, 0.75). In that case the hnear statis- 
tics performed poorly and, as one would expect, were dominated by the quadratic 
statistics. 



5 Simulation Studies for Case-Control Settings 

Here, we provide detailed numerical results for case-control studies when a normal 
approximation for S might not be adequate. As in previous sections, we examine the 
performance of W^s^ ^lw^ Wc and Wh- We first considered a normal approximation 



for the linear statistics and linear combination of chi-squares (see (2.7)) for quadratic 
statistics for obtaining p-values or Type I errors, and investigated their adequacy 
by simulation. We then conducted simulations to assess the statistics power under 
different scenarios. 

We assume that the distribution of Yi given Xi in the population is Bernoulli with 

and that the Xij in the population are mutually independent Bernoulli variables 
with P{Xij = 1) = pj for j = 1, J. Similar to Section 4, we consider for power 
assessments a broad range of scenarios by randomly generating 500 different models. 
Each model has randomly selected values J, Jc, Jd, Pj and pj as follows: J was 
randomly selected from the set {10,20,30,40,50}; pc = Jc/J was generated from 
f/(0.1,l); pd = Jd/Jc was generated from ?7(0.75,1), /3j was generated differently 
depending on the scenarios described later and pj was generated from f/(0.005, 0.02). 
Sample size for both cases and controls was fixed to be ni = uq = 500, the level of the 
test was set to be a = 0.0001 and without loss of generality we took /3o = —2.1922 
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(giving P{Yi = l\Xi = 0) = 0.1). We consider two different situations concerning pj 
and Pj, which is a log odds ratio, log(OR). The first scenario assumes no relationship 
between pj and (3j for each SNP and takes e'^^' ~ f/(1.5,3) for each causal SNP. The 
second scenario assumes that the odds ratio for a causal SNP is inversely proportional 
to ^/pj{l - Pj); we set el'^^l = C / y^pj{l - Pj), where C = 4^0.005(1 - 0.005). Under 
these scenarios the odds ratio for deleterious variants ranges from 2 for pj = 0.02 to 
4 for Pj = 0.005 and from 0.5 down to 0.25 for protective variants. 

Due to the sparsity of rare variants and the dichotomous trait, normal and chi- 
square approximations for statistics may produce inflated or conservative Type I 



errors as was observed in previous studies (Lin and Tang (2011), Basu and Pan 



(2011)). To investigate the suitability of the approximations for nominal Type I 
error a = 0.05,0.01,0.001 and 0.0001, we considered the model with pj = 0.01 for 
all J SNPs under null hypothesis (/3 = 0). Results are displayed in Table 1 for 
linear statistic Wis and Hotelling statistic Wh (Table 1). Empirical Type I errors in 
the table are based on 10^ simulation replicates, and are the proportions of the 10^ 
samples for which the corresponding test statistic exceeded the a critical values given 
by a normal or chi-square approximation. The 10^ simulation replicates could be 
realized efficiently because the Sj are mutually independent and are functions of the 
Tj and rrij only; under Ho, rrij is generated from Binomial{1000, 0.01) and Tj from 



(3.2). Table 1 shows that the normal approximation for Wl is accurate, but quadratic 
tests have conservative Type I errors with chi-square approximations. Therefore, for 
each randomly chosen scenario for assessing power, we first generated 10^ models 
under the null hypothesis (/3 = 0), to obtain empirical critical values for Type I error 
a = 0.0001. 

The case-control samples under Hi were simulated using for computational conve- 
nience the assumption that the Xij are mutually independent, given = or yj = 1. 



23 



Under this assumption, Xi = {Xn, ...,Xijy, for the i^^ individual with case status is 
generated from Xij = Binomial{l, p{Xij = l\Yi = 1)). Similarly, X^ = {Xn, ...,Xij)' 
for i^^ individual with control status is generated from Xij = Binomial{l^p{Xij = 
l\Yi = 0)). Conditional probabilities of Xij = 1 given = or 1^ = 1 were calculated 
using Monte-Carlo sampling. This is an approximation, but empirical powers were 
found to be very close to those obtained under much more computationally intensive 
simulations that are based on the exact model. For each combination of parameters, 
simulations to assess power consisted of 1,000 replications, with power estimated by 
the proportion of samples for which test statistic exceeded its critical value. 

Figure 7 shows results based on the 500 randomly generated models under the first 
scenario. It compares the empirically calculated power for linear statistics with and 
without weighting, Wls and Wlw- It also compares power of the Hotelling statistic 
Wh and C-alpha statistic Wc- We see that the C-alpha performs better than the 
Hotelling statistic, while the linear statistic with no weights has an advantage over 
the statistic with weights. Under case-control scenarios, deleterious SNPs have larger 

n 

values of = ^ Xij than neutral SNPs with the same pj due to enrichment from 

i=l 

the cases. This explains why Wls and Wc have an advantage over Wlw and Wh, 



respectively. Similar results were obtained by Basu and Pan (2011). 

Figure 8 shows results based on the 500 randomly generated models under the 
second scenario. We again see a difference in power between linear statistics Wls 
and Wlw- However, the systematic power difference between the quadratic statistics 



is absent in Figure 8. This supplements the results of Basu and Pan (2011), who 
did not consider cases where the genetic effect is inversely proportional to MAF, and 
indicates that the relative performance of Wc and Wr depends on the relationship 
between SNP effects and MAF. 

We also compare power of the linear statistic Wls and quadratic statistic Wc for 
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the 500 models under the first scenario in Figure 9, and Figure 10 compares power 
of those statistics under the second scenario. Similar conclusions to those for the 
quantative trait study can be made. The variation in power between the two tests is 
quite large here but in Figure 10 especially, the linear statistics achieves high power 
more often. We also investigated the relationships between model parameters and 
the power of the linear (Wis) and quadratic (Wc) statistics. We present results for 
settings with J = 30 and rii = uq = 500 in Figure 11 under assumption of independent 
of Pj and Pj and in Figure 12 under the odds ratio for a causal SNP is inversely 
proportional to ^/pj{l — Pj)- Similarly to numerical comparisons under normality, 
we observed that linear statistics have lower power than quadratic statistics even 
when all or almost all causal SNPs have effect in the same direction. One difference 
from results in Section 4 is that linear statistics has smaller drop of power when large 
proportion of causal SNPs have effect in the different direction. This is because the 
association statistics Sj in the case-control setting is scaled MAF difference between 
cases and controls. Under the logistic model, this difference is not symmetric for 
deleterious and protective SNPs even when they have same effect but different sign. 



6 Applications: GAW17 Data 



The numerical studies in Sections 4 and 5 assume that SNPs were mutually inde- 
pendent. To consider settings were this might no be so, we examined real human 



sequence data (1000 Genomes Project Consortium (2010)) thet were used to generate 
quantitative trait data in Genetic Analysis Workshop 17(GAW 17); see Almsay et al. 
(2011). 

We analyzed quantitative trait Q2 which was influenced by 72 SNPs (70 with MAF 
in the range (0.16%, 17%)) in 13 genes (see Table 2 of Almasy et al. ( |2011 )) but not 
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by other covariates . We used data from the n = 321 unrelated Asian subjects (Han 
Chinese, Denver Chinese and Japanese). We calculated permutation-based p- values 
for the four tests, Wis, Wlw, Wc and Wh discussed in previous sections, and we 
estimated power by analyzing all 200 replicates provided by GAW17 at the a = 0.05 
level (Table 2). The choice a = 0.05 was based on the overall low power for detecting 
effects due to small sample size, small genetic effect, extremely small MAF and the 
low proportion of causal variants in a gene. 

Results in Table 2 are consistent with our previous conclusions: (i) linear tests 
with and without weighting based on MAF vary in relative power; (ii) quadratic 
statistics Wc and Wh also have variable relative power; (iii) relative performance of 
the linear and quadratic tests is highly variable. As expected, with a large number of 
causal SNPs linear statistics can outperform quadratic (e.g. genes SIRTl, SREBFl), 
but when the proportion of causal SNPs is low quadratic statistics outperform linear 
(e.g. genes BCHE, PDGFD, RARE). 



7 Regression Models and Additional Covariates 



In some settings a test of no association may be based on a regression model (e.g. 



Morris and Zeggini] (|2010j)). In fact, |Lin and Tang| (|2011|) and |Wu et al.| ( |2011[ ) have 
stressed that adjustment for covariates and population stratification will be important 
in many contexts involving rare variants. We now discuss how our framework is 
extended to deal with the inclusion of covariates; for illustration, we consider the case 
of a binary trait. Suppose that in addition to the genotype vector Xi there is a vector 
Vi of covariates that may be related to a binary trait Yi. Then a logistic regression 
model 

exp(/3o + /3'X, + 7V,) 



Pr(y, = 1|X„*;,) = 



1 + exp(/3o + f3'X, + j'v, 
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(7.i: 



might be considered, and a test of Hq : f3 = can be carried out. For testing rare 

J 

variants it is common to replace the term f3'Xi in (7.1) with /3rj, where rj = ^ Xij 



is the total number of rare variants (e.g. Morris and Zeggini (2010); Yilmaz and Bull 



(2011)), but this corresponds to using a linear statistic in previous sections and is 



often ineffective. We consider the case where f3 = . . . , in order to examine 
settings for which causal SNPs may be either deleterious or beneficial. In that case 
consideration of the power of alternative tests in large samples parallels the discussion 
in Section 4, as follows. 

Let f3 be the estimator of f3 based on the model in question and assume that under 
Hq : f3 = 0, the asymptotic distribution of y/n^ is multivariate normal with mean 



and covariate matrix Following Li and Lagakos (2006), we consider a sequence of 
contiguous alternatives 

i/J") : f3 = hj^ (7.2) 

where h = (6i,...,6j)' is a specified vector. Under this sequence as n — )■ oo the 
distribution of ^/nf3 approaches a multivariate normal distribution with mean b and 
covariance matrix Thus, asymptotic power for a test statistic can be computed 



in the same way as in Section 4. Li and Lagakos (2006) compare the "omnibus" test 



— 1 



statistic W = (3^ (3, where X] is a consistent estimate of ^ under Hq, with linear 



statistics Z = a'f3. These are analogous to (2.5) and (2.4), respectively. In fact 



note that if we consider the linear regression model (4.1) with the independent 



iV(0, a^) random variables, then ^ = {X'X)-^X'Y, where 1^ = (^i, . . . , F„)' and X 
is a centered n x J matrix whose i'th row is Xi = {Xn, Xjj)' — {Xi, . . . , Xj), and 
so 



-1 



where S has j'th element Sj = ^ (Xjj — Xj) Yi = Xij (Yi — F) . 



1=1 



i=l 
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A similar results hold in the logistic model (7.1). The "omnibus" Wald statistic 



W* = ^Yl ^ is asymptotically equivalent to the likelihood score statistic Ws for 
testing /3 = and when there are no covariates Vi, Ws is exactly Wh = S'T.g'^S. 
For the case where covariates Vi are present, an asymptotically equivalent statistic to 



W* is based on the likelihood score for prospective sampling under (7.1). The score 
statistic for testing /3 = is easily found as 



(7.3) 



1=1 



where jdi = e^°^^ / (1 + e^°^'^'^') and (3o, 7 are estimated from (7.1) when 



/3 = 0. It also follows from standard maximum likelihood large sample theory that 
the covariance matrix of U under Hq is estimated consistently by 



i=l 



.i=l 



i=l 



(7.4) 



where af = fii{l — jli) and Vi = {l,v'^'. These correspond to results given by 



Lin and Tang (2011), who consider linear statistics based on linear combinations 



of the elements Ui,...,Uj of U. The statistic (7.3) and variance estimate (7.4) can 



be shown to apply under case-control sampling and they give test statistics such as 
W;j = U'Jl^^U and Wl = {w'U)/{w'i:^^w), which correspond to Wh and Wl in 



preceding sections. When there are no covariates Vi, it is readily seen that (7.3) 



reduces to (3.1) and that (7.4) equals (n — l)/n times (3.3). It should be noted that 



when covariates Vi are present, the normal approximations considered earlier apply, 
but the permutation distribution p-values do not unless the Xi are independent of 



the Vi. Lin and Tang (2011) suggest a parametric bootstrap as an alternative, based 



on randomly generating responses from the fitted null model based on Pq, 7. We note 
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this can be computationally intensive in case-control settings 

It should be mentioned that in the case of quantitative Y-dependent sampling 



and models with supplementary covariates Vi as in (7.1), adjustments to estimating 
functions (e.g. based on inverse probability of selection weights) are needed; this is 
beyond our present scope. 

8 Discussion 

In this paper, we have compared tests of association between rare variants and phe- 
notypes within a unified framework which gives theoretical insights about the perfor- 
mance of the methods. Our treatment handles all types of phenotypes and can deal 
with phenotype-dependent sampling of individuals and correlation among the SNPs 
under consideration. One of the important conclusions of this work is that depending 
on the alternative, methods can have greatly varying power. We found in extensive 
numerical studies that, as expected, when both deleterious and protective SNPs are 
present the quadratic test statistics are much better. They also outperform linear 
statistics in settings where causal SNPs are all deleterious (or all protective), but a 
substantial fraction of the SNPs are neutral (not associated with the phenotype). Al- 
though they perform well across a broader range of settings, relative performance can 
vary according to the fraction of causal SNPs and importantly, according to whether 
causal affects are associated with lower MAF. Our results also indicate that power to 
detect moderate levels of association is not high unless sample sizes are very large or 
a high proportion of the chosen SNPs are causal. Consequently it is critical to ob- 
tain relevant biological information that can guide the selection of SNPs or weighting 



strategy for pooled association testing (e.g. King et al. (2010)). 



Our results supplement these of Basu and Pan (2011), and a brief comparison is 
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useful. They found similar results to ours in simulation studies for case-control scenar- 
ios, concerning the performance of linear statistics. Among the quadratic statistics, 
they found that C-alpha/SSU type statistics Wc = S'S was generally best, and su- 



perior to the Hotelling statistics (2.5). However, their simulation scenarios did not 
include cases where causal effects were associated with SNPs having smaller MAFs. 
Our numerical studies and investigation of GAW17 data indicate the importance of 
relationships between SNP effects and MAFs. As an approach to rare variant test- 
ing in the absence of strong prior information, we support the recommendation of 



Basu and Pan (2011), to perform tests using both linear and quadratic statistics. In 



Derkach et al. (2012), we examine statistics that combine evidence from linear and 



quadratic tests. 
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Figures and Tables 

Comparison of linear statistics Comparison of quadratic statistics 
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Figure 1: 

Comparison of power of linear statistics: line ar st atistic Wis is given by (2.4) with 



set to 1 and linear statistic Wlw given by (2.4) with wj set to 1/ vPj(l — Pj 



Comparison of power of quadratic statistics: Hoteiiing statistic is Wh (2.5) and C- 
alpha statistic is Wq (2.6) with matrix A set to the identity /. 

The points in each plot show the power of each statistic for each of the 1000 randomly 
chosen models with the genetic effect of a SNP independent of its MAF. 
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Comparison of linear statistics Comparison of quadratic statistics 
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Figure 2: 

Comparison of power of linear statistics: line ar st atistic Wls is given by (2.4) with 



Wj set to 1 and linear statistic Wlw given by (2.4) with Wj set to —Pj)- 

Comparison of power of quadratic statistics: Hoteiiing statistic is Wh (2.5) and C- 
alpha statistic is Wq (2.6) with matrix A set to the identity /. 

The points in each plot show the power of each statistic for each of the 1000 randomly 
chosen models with with explained variation of a single SNP independent of its MAF. 
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Figure 3: 

Comparison of power of linear and quadratic statistics: Wl (2.4) with weights Wj set 
to 1 (Wls) and the Hotelhng statistic Wh ( [2^ 

Points represent the power for each of the 1000 randomly chosen models with param- 
eters described in Section 4 and with genetic effect of a SNP independent of its MAF. 
n - sample size. 
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Figure 4: 

Comparison of power of linear and quadratic statistics: Wl (2.4) with weights Wj set 
to 1 (Wis) and the Hotelling statistic Wh ( [2^ 

Points represent the power for each of the 1000 randomly chosen models with param- 
eters described in Section 4 and with explained variation of a single SNP independent 
of its MAF. n - sample size. 
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(b) 



Figure 5: Empirical power of the statistics represented as function of two parameters: 
J-number of SNPs, n = 1000 

On the X-axis: number of causal SNPs, larger scale; number of deleterious SNPs, 
smaller scale 

(a) Linear statistic, (b) Hotclling statistic. 

Results are based on 10000 randomly chosen models with parameters described in 
Section 4, with the genetic effects of SNPs independent of MAP. 
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X-axes are 
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Smaller scale: number of deleterious 
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(b) 



Figure 6: Empirical power of the statistics represented as function of two parameters: 
J-number of SNPs, n = 1000 

On the X-axis: number of causal SNPs, larger scale; number of deleterious SNPs, 
smaller scale 

(a) Linear statistic, (b) Hotclling statistic. 

Results are based on 10000 randomly chosen models with parameters described in 
Section 4, with explained variation of single SNPs independent of MAP. 



38 



a ^ 0.05 



a = 0.01 



Methods 



Methods 



J 


LS 




QS(Chsq) 


LS QS(Chsq) 


10 


0.049 




0.041 


0.0098 0.006 


20 


0.050 




0.041 


0.0099 0.006 


50 


0.050 




0.041 


0.0099 0.007 


100 


0.050 




0.041 


0.0100 0.007 




a = l 


10 


-3 


a = l- 10"^ 




Methods 


Methods 


J 


LS 




QS(Chsq) 


LF QS(Chsq) 


10 


0.94 • 10- 


-3 


0.3 • 10-3 


0.84- 10-* 0.14-10-4 


20 


0.99 • 10- 


-3 


0.4- 10-3 


0.88-10-4 0.15-10-4 


50 


1.01 • 10- 


-3 


0.4- 10-3 


1.02 - 10-4 0.38 - 10-4 


100 


1.01 • 10- 


-3 


0.5 • 10-3 


0.99 - 10-4 0.48 - 10-4 



Table 1: Empirical type 1 errors for two tests of association with a binary outcome. LS: normal 



quadratic statistic (2.5) 



approximation for the linear statistic (2.4|; QS(Chsq): chi-square approximation for the Hotelling 



Sample size is 500 cases and 500 controls; MAF of each SNP is 0.5%; the number of simulation 
replicates is 10^. 
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Linear statistics 



Quadratic statistics 
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Figure 7: 

Comparison of power between linear statistics: (2.4) with weights Wj set to 1 (Wls) 
and with weights Wj set to 1/ ^/pj{l — Pj) (Wlw)- 



Comparison of power between quadratic statistics: Hotelhng statistic Wh (2.5) and 
C-alpha statistic Wq (2.6) with matrix A set to the identity /. 

Points represent the power for each of the 1000 randomly chosen models with param- 
eters as described in Section 5, with the genetic effect of a SNP f3j independent of its 
MAF pj. 
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Linear statistics Quadratic statistics 




Power of linear statistic withiout weights Power of Hotelling statistic 



Figure 8: 

Comparison of power between linear statistics: (2.4) with weights Wj set to 1 (Wls) 
and with weights Wj set to 1/ ^ypj{l — Pj) (Wlw)- 

Comparison of power between quadratic statistics: Hotelhng statistic Wh (2.5) and 
C-alpha statistic Wq (2.6) with matrix A set to the identity /. 

Points represent the power for each of the 1000 randomly chosen models with param- 
eters as described in Section 5, with el'^-'l = C/ ^/pj{l — Pj)- 
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Linear vs Quadratic statistic 
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Figure 9: 



Comparison of power between linear statistic (2.4) with weights Wj set to 1 (Wls 



and C-alpha statistic Wq (2.6) with matrix A set to the identity /. 

Points represent the power for each of the 1000 randomly chosen models with param- 
eters as described in Section 5, with the genetic effect of a SNP /3j independent of its 
MAF pj. 
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Linear vs Quadratic statistic 



« • • • •• 



. • • • 



0.0 



0.2 0.4 0.6 

Power of quadratic statistic 



0.8 



1.0 



Figure 10: 



Comparison of power between linear statistic (2.4) with weights Wj set to 1 (Wls 



and C-alpha statistic Wq (2.6) with matrix A set to the identity /. 

Points represent the power for each of the 1000 randomly chosen models with param- 
eters as described in Section 5, with parameters described in Section 5 and with with 
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(b) 



Figure 11: Empirical power of the statistics represented as function of two parame- 
ters: J-number of SNPs, n = 1000 

On the X-axis: number of causal SNPs, larger scale; number of deleterious SNPs, 
smaller scale 

(a) Linear statistic, (b) HotcUing statistic. 

Results are based on 500 randomly chosen models with parameters described in Sec- 
tion 5, with the genetic effect of a SNP /3j independent of its MAF pj . 
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X-axes are 
Larger scale: number of causal 
• Smaller scale: number of deleterious 
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(b) 



Figure 12: Empirical power of the statistics represented as function of two parame- 
ters: J-number of SNPs, n = 1000 

On the X-axis: number of causal SNPs, larger scale; number of deleterious SNPs, 
smaller scale 

(a) Linear statistic, (b) HotcUing statistic. 

Results are based on 500 randomly chosen models with parameters described in Sec- 
tion 5 and with with e^^^^ — C/ \/Pj{^ ~ Pj)- 
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Supplementary Materials 



Derivations 

Proof of (2.3) with discrete Y 

In this case, let Y{, . . . , denote the K distinct values of 1^ {i = 1, . . . ,n), and 

k 

let n* be the number of values Y*, with — n. Let a* be a rank score assigned 

r=l 

k 

to Y*{r = 1, k), centered so that ^ n*a* = 0. Then under Hq, we have 

r=l 



Pr (a, = a*) = — , 
n 

Pr {ai = a^/ = a*) = ^'^ , ~ , i ^ i' 

n[n — 1) 

Thus E{ai) — and for i ^ i', 



k k 



E{aiai') = q:*q:* Pr (aj = a*, aji = a*) 



^ n(n-l) A^A^ ^ ^n(n-l) 

r=l ^ ' r^s 

~ ^^nin-l) ^ n(n - 1) 

r=l s=l ^ ' r=\ ^ ' 
n 9 

^ n(n — 1) 

1=1 ^ ' 
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Thus Ep{Sj) = and 

n n 



i=i i'=i 

n 



i=l i^i' 
n / n 2\ /" 2 



which reduces to (2.3). 
Proof of (3.3) 

Since Ep{Sj) — the covariance of Sj and 5'^ is 

n n 

Ep(SjSe) - E E (^^^^') i^'^ - i^^'i - ' 
i=i i'=i 



using the fact that (3.1) also equals ^{Xij — Xj)Yi. Thus 



Ep{SjSe) = J2Ep{Y^){X,j-Xj){X,e-Xe) 

1=1 

+ E E^^ (^^^^') i^i^ - i^i'i - 



n 
i=l 

^ ^ i=l 



Noting that ^Xij — ruj and that XijXig^ — rriji, we obtain (3.5); noting that 

i=l i=l 

n 

^ Xfj = rrij, we also obtain (3.3). 



i=l 
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